In [7]:
f_esfera = lambda xk: xk[0]**2 + xk[1]**2
f_esfera_grad = nd.Gradient(f_esfera)
f_esfera_hessian = nd.Hessian(f_esfera)
plot_function(f_esfera, title = 'Função Esfera', ndim = 3)
# importando funçõs necessárias
!pip install numdifftools
from IPython.display import clear_output
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from matplotlib import cm
import numdifftools as nd
import warnings
clear_output()
pd.set_option('display.float_format', lambda x: '%.8f' % x)
np.set_printoptions(8)
warnings.filterwarnings("ignore")
# funções auxiliares
def plot_function(f, title, ndim, dom = np.linspace(-5, 5, 500), angle = (20, 20)):
"""
Função que plota gráficos em 2-d ou 3-d.
Parâmetros
----------------------------------------
f : callable
Função objetivo (função custo).
title : str
Título do gráfico
ndim : str (2 ou 3)
Dimensão da função
dom : list or np.arange or np.linspace, opcional
Domínio da função desejado
angle : tuple, opcional
Ângulo para ver o gráfico. Padrão = (20, 20)
Saída
----------------------------------------
plot
"""
plt.style.use('fivethirtyeight')
# 2-d plot: y = f(x)
if ndim == 2:
X = dom
Y = f(X)
fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='2d')
ax.set_title(title)
ax.set_xlabel('$x$')
ax.set_ylabel('$f(x)$')
ax.plot(X, Y)
ax.view_init(angle[0], angle[1])
plt.tight_layout()
plt.show()
#3-d plot: z = f(x, y)
else:
x, y = dom, dom
X, Y = np.meshgrid(x, y)
Z = f([X, Y])
fig = plt.figure(figsize=(8,8))
ax = fig.gca(projection='3d')
ax.set_title(title)
ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$f(x_1, x_2)$')
ax.plot_surface(X, Y, Z, cmap='jet')
ax.view_init(angle[0], angle[1])
plt.tight_layout()
plt.show()
def plot_results(f, xi, yi, method, dim = [np.linspace(-7, 7, 500), np.linspace(-7, 7, 500)]):
"""
Função que plota os result
Parâmetros
----------------------------------------
f : callable
Função objetivo (função custo).
xi : np.array
Valores de x em cada iteração do método
yi : np.array
Valores de y em cada iteração do método
method: str
Método desejado: '1', '2' ou '3'
dim : np.linspace, opcional
Valores de x e y para para plotar as curvas de nível
Saída
----------------------------------------
plot
"""
plt.style.use('fivethirtyeight')
if method.lower().strip() == '1': title = 'Descida do gradiente com busca de Armijo'
elif method.lower().strip() == '2': title = 'Método de Newton com busca de Armijo'
else: title = 'Método BFGS com busca de Armijo'
# z = f(x,y)
if xi.shape[1] == 2:
x, y = dim[0], dim[1]
X, Y = np.meshgrid(x, y)
Z = f([X, Y])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 8))
plt.suptitle(title, y = 1.05)
# curvas de nível
ax1.plot(xi[:,0], xi[:,1], linestyle='--', marker='o', color='black', linewidth = 3)
ax1.plot(xi[-1,0], xi[-1,1], 'ro', markersize = 11)
ax1.set(title='Caminho durante a otimização - Curvas de Nível', xlabel='x1', ylabel='x2')
CS = ax1.contour(X, Y, Z, 15, cmap = 'jet')
ax1.clabel(CS, fontsize='smaller', fmt='%1.2f')
# valor da função custo em cada iteração
ax2.plot(yi, linestyle='--', marker='o', color='black')
ax2.plot(len(yi)-1, yi[-1], 'ro', markersize = 11)
ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
ax2.set(title = 'Valor da função custo durante a otimização', xlabel='Iterações', ylabel='Valor da função custo')
ax2.legend(['Busca de Armijo'])
plt.show()
# y = f(x)
else:
x = dim[0]
y = f(x)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))
plt.suptitle(title, y = 1.05)
ax1.plot(xi[:,0], xi[:,1], linestyle='--', marker='o', color='black', linewidth = 3)
ax1.plot(xi[-1,0], xi[-1,1], 'ro', markersize = 11)
ax1.set(title='Caminho durante a otimização', xlabel='x1', ylabel='x2')
ax2.plot(yi, linestyle='--', marker='o', color='black')
ax2.plot(len(yi)-1, yi[-1], 'ro', markersize = 11)
ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
ax2.set(title = 'Valor da função custo durante a otimização', xlabel='Iterações', ylabel='Valor da função custo')
ax2.legend(['Busca de Armijo'])
plt.tight_layout()
plt.show()
def printa_resultados(chute, xi, yi, grad_norm):
chute = np.round(chute, 8)
print('\n*************** Resultados ***************')
print()
print('Número de iterações: ', len(yi))
print()
print('***** Primeira iteração *****')
print()
print(f'f({chute}) = {np.round(yi[0], 8)}')
print()
print(f"||∇f({chute})|| = ", np.round(grad_norm[0], 8))
print()
print()
print(f'***** Última iteração: {len(yi)}° *****')
print()
print(f'f({np.round(xi[-1], 8)}) = {np.round(yi[-1], 8)}')
print()
print(f"||∇f({np.round(xi[-1], 8)})|| = ", np.round(grad_norm[-1], 8))
def minimiza_f(f, metodo, chute, sigma = 0.02, tol = 1e-8):
""" Função para testar os métodos
Parâmetros
----------------------------------------
f : callable
função objetivo no formato estabelecido
metodo : int (1, 2 ou 3)
1 - Método do Gradiente
2 - Método de Newton
3 - Método BFGS
chute : np.array ou list
array com o ponto inicial (mesma dimensão da função)
sigma : float > 0, opcional
cte de decréscimo de Armijo. Padrão - 0.02
tol : float > 0, opcional
tolerância pré-definida. Padrão = 1e-8
Saída
----------------------------------------
Resultados obtidos da minimização """
f_grad = nd.Gradient(f)
m = {1:'Método de Descida do Gradiente com condição de Armijo', 2:'Método de Newton com condição de Armijo', 3: 'Método Broyden–Fletcher–Goldfarb–Shanno (BFGS) com condição de Armijo'}
print(f"{m[metodo]} \n\n")
if metodo == 1:
xi, yi, grad_norm = GradientDescent(f, f_grad, chute = chute, sigma = sigma, tol=tol)
printa_resultados(chute, xi, yi, grad_norm)
elif metodo == 2:
f_hessian = nd.Hessian(f)
xi, yi, grad_norm = Newton(f, f_grad, f_hessian, chute = chute, sigma = sigma, tol=tol)
if len(xi) != 0:
printa_resultados(chute, xi, yi, grad_norm)
else:
xi, yi, grad_norm = BFGS(f, f_grad, chute = chute, sigma = sigma, tol=tol)
printa_resultados(chute, xi, yi, grad_norm)
if len(chute) == 2 and len(xi) != 0:
print('\n\n\n\n')
plot_results(f, xi, yi, dim = [np.linspace(int(np.min(xi[:,0]))-7, int(np.max(xi[:,0]))+7, 1000), np.linspace(int(np.min(xi[:,1]))-7, int(np.max(xi[:,1]))+7, 1000)], method = str(metodo))
Uma das principais aplicações do Cálculo Diferencial e Integral é na área de otimização, uma vez que, na maioria das vezes, o objetivo é maximizar ou minimizar determinada função.
Alguns exemplos disso vão desde a otimização de carteira, no mercado financeiro; caminhos, para serviços de entrega e sistemas GPS; e em diversos modelos de aprendizado de máquina, nos quais a função custo deve ser minimizada para se obter os melhores parâmetros que representarão adequadamente seus dados para determinada tarefa.
Nisso, a otimização não linear apresenta várias maneiras de lidar com tais problemas. Neste trabalho serão apresentados os métodos clássicos de descida, sendo eles: Método do gradiente, de Newton e BFGS, todos com busca linear de Armijo. Junto com cada método será apresentada sua descrição teórica, de algoritmo e uma breve análise da minimização de algumas funções, que foram retiradas de exemplos dados em $[1]$, $[2]$ e $[3]$. A parte teórica foi inspirada em notas de aula, em $[4]$ e $[5]$.
A linguagem de programação escolhida foi Python, versão 3.7.13.
!python --version
Python 3.9.7
A Condição de Armijo (1° Condição de Wolfe) é definida da seguinte maneira:
$$ $$$$ \begin{equation}\tag{1} f(x_{k} + \lambda_{k}d_{k}) < f(x_{k}) + \sigma\lambda_{k}\nabla^{T}{f(x_k)}d_{k} \end{equation} $$$$ $$Onde:
$x_k \equiv \text{chute inicial (initial guess)}$
$\sigma \in (0, 1) \equiv \text{constante de decréscimo}$
$\lambda_{k} > 0 \equiv \text{passo na iteração } k$
$\nabla^{T}{f(x_{k})}\,d_{k} < 0 \equiv \text{direção de descida em } x_{k}$
$$ $$Assim, a condição estabelece uma maneira de calcular o passo necessário para um Algoritmo de Descida, como os que serão apresentados abaixo. Além disso, vale ressaltar que este decréscimo da função é proporcional ao passo. Dessa maneira, esquematizando o Algoritmo de busca com condição de Armijo, temos:
$$ $$Algoritmo 1 - Busca Linear com Condição de Armijo
Dada uma direção de descida $d_{k}$:
$$ $$Importante ressaltar que $\gamma$ é uma constante em que $\gamma \in (0, 1)$.
Abaixo está a implementação do Algoritmo de busca com condição de Armijo (Algoritmo 1).
def Armijo_Search(f, xk, dk, sigma = 0.02, gamma = 0.5):
'''
Busca linear com Condição de Armijo (1° Condição de Wolfe) - Desigualdade 1, que consiste em fazer um decréscimo da f proporcional ao tamanho do passo
Parâmetros
----------------------------------------
f : callable
Função objetivo (função custo).
xk : array
Ponto atual.
dk : array
Direção de descida (grad_T(f(x_{k})) . d_{k} < 0).
sigma : float, opcional
Valor de sigma entre (0, 1) - constante de decréscimo. Padrão = 0.02
gamma : float, opcional
Valor de gamma entre (0, 1). Padrão = 0.5
Saída
----------------------------------------
lambda : float
Valor de lambda que satisfaz a condição de Armijo.
f_x0 : float
Valor de f no ponto x_{k+1}.
'''
# Começa constante lambda = 1
lambda_ = 1
# Direção de descida d = grad_{T}(f(x_{k})) . dk < 0
desc_direction = np.dot(nd.Gradient(f)(xk), dk)
f_x0 = f(xk + lambda_ * dk)
# Calcula a função em um ponto menor que x_k com passo inicial lambda = 1
# Condição de Armijo para determinar o tamanho do passo, diminuindo a função custo
while (f_x0 > f(xk) + sigma * lambda_ * desc_direction):
lambda_ = gamma * lambda_
f_x0 = f(xk + lambda_ * dk)
# Retorna o passo lambda e o valor da função
return lambda_, f_x0
O Método do Gradiente pode ser de dois tipos: Subida ou Descida. O interesse nesse trabalho é o de minimização de funções, logo, o de Descida do Gradiente, uma vez que o vetor gradiente fornece a direção de maior crescimento da função.
Esse método consiste em se obter uma direção $d_k$ tal que $d_k = - \nabla f(\vec{x})$ e a partir disso calcular o novo ponto $x_{k+1}$ em que a norma do gradiente $\lVert \nabla f(x_{k+1}) \rVert$ seja menor e o custo da função $f(x_{k+1})$ diminua. O algoritmo esquematizado está abaixo:
Algoritmo 2 - Método do Gradiente com condição de Armijo
Dado um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$) e para $\nabla f(x_{k}) \neq 0$, realizar os seguintes passos:
1. Calcular a direção de descida $d_{k} = - \nabla f(x_{0})$;
2. Realizar Busca de Armijo e encontrar um $\lambda_k$ que satisfaça a desigualdade 1;
3. Calcular o novo ponto $x_{k+1} = x_{k} + \lambda_{k} d_{k} = x_{k} - \lambda_{k} \nabla f(x_{k})$;
Em que tais etapas são realizadas até algum critério de parada pré-estabelecido. No presente trabalho foram 2 critérios: ou a parada se dava por um alguma tolerância referente à norma do gradiente da função ($\rVert\nabla f(x_{k}) \lVert \leq 10^{-18}$), por exemplo, ou por um número máximo de iterações (500 por $\textit{default}$). Ambos os parâmetros são alteráveis nas funções apresentadas.
A convergência do Método de Descida do Gradiente é linear, ou seja, a sequência {$x_k$} converge para $x^{*}$ quando:
\begin{equation} \lVert x_{k+1} - x^{*} \rVert \leq \epsilon \, \lVert x_{k} - x^{*} \rVert \end{equation}$\forall \, k \in \mathbb{N}$ e com $\epsilon > 0$.
A chamada para a Função que implementa o Método do Gradiente com Condição de Armijo é:
GradientDescent(f, f_grad, chute, sigma, tol)
Onde:
$$ $$f é a função desejada (tipo: callable)
f_grad é o gradiente da função f (tipo: callable)
chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)
sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)
tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = $10^{-8}$)
E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.
$$ $$x são pontos $x_k$ calculados pelo Método
y são valores de $f(x_{k})$ calculados pelo Método
grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método
Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:
$$ $$x, y, _ = GradientDescent(f, f_grad, chute, sigma, tol)
Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 2:
def GradientDescent(f, f_grad, chute, sigma = 0.02, tol=1e-8):
"""Algoritmo de Descida do Gradiente com busca linear com condição de Armijo. Esquema:
I - Direção de descida: d_k := −∇f(x).
II - Determinação do passo com busca de Armijo.
III - Obtém o próximo candidato.
Parâmetros
--------------------
f : callable
Função custo
f_grad : callable
Gradiente da função f
chute : array
Valor inicial de x ("chute")
sigma : float, opcional
Constante de Decréscimo de Armijo. Padrão = 0.02
tol : float, opcional
Tolerância padrão de 1e-8
Saída
--------------------
xk : array
valores de xk
yk : array
valores da função custo em cada xk
grad_f : array
valores da norma do gradiente da f em xk """
# Valores iniciais de xk, fk e grad_fk
xk = chute
fk = f(xk)
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
max_iter = 500
# Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
num_iter = 0
x_pontos = [xk]
y_pontos = [fk]
grad_f = [grad_fk_norm]
print(f'Chute inicial: y = {fk}, x = {xk} \n')
# Calcula nova iteração com busca de Armijo
while (grad_fk_norm > tol and num_iter < max_iter):
# Determina a direção
pk = -grad_fk
if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16):
print("\nErro: Overflow\n")
break
# Faz a busca de Armijo, obtendo o passo lambda e a função custo naquele passo
lambda_, fk = Armijo_Search(f, xk, pk, sigma)
# calcula x_{k+1}
xk = xk + lambda_ * pk
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
# Itera mais uma vez
num_iter += 1
x_pontos.append(xk)
y_pontos.append(fk)
grad_f.append(grad_fk_norm)
# print results
if num_iter == max_iter:
print('\nNúmero de iterações máximo atingido.\n')
return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)
Nesta parte algumas funções exemplo serão testadas utilizando o Método de Descida do Gradiente com busca linear de Armijo.
O primeiro exemplo mais simples é com a Função Esfera, definida como:
$$ $$\begin{equation} f(\mathbf{x}) = \sum_{i=1}^{n} x^{2}_{i} \end{equation}$$ $$Com $\mathbf{x} \in \mathbb{R}^{n}$, $n \in \mathbb{N}$. Foi considerado para os exemplos abaixo a função com 2 variáveis, ou seja, $n=2$.
f_esfera = lambda xk: xk[0]**2 + xk[1]**2
f_esfera_grad = nd.Gradient(f_esfera)
f_esfera_hessian = nd.Hessian(f_esfera)
plot_function(f_esfera, title = 'Função Esfera', ndim = 3)
Cujos resultados são os seguintes:
chute = [4, -3]
xi, yi, grad_norm = GradientDescent(f_esfera, f_esfera_grad, chute = chute, tol=1e-10)
df1_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
df1_1
Chute inicial: y = 25, x = [4, -3]
| x1 | x2 | f(x1,x2) | ||∇f|| | |
|---|---|---|---|---|
| 0 | 4.00000000 | -3.00000000 | 25.00000000 | 10.00000000 |
| 1 | 0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 |
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 2 ***** Primeira iteração ***** f([ 4 -3]) = 25.0 ||∇f([ 4 -3])|| = 10.0 ***** Última iteração: 2° ***** f([0. 0.]) = 0.0 ||∇f([0. 0.])|| = 0.0
O algoritmo convergiu de acordo com critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$, estabelecido anteriormente com $2$ iterações.
Abaixo temos 2 $\textit{plots}$, o primeiro com as curvas de nível e o caminho tomado pelo algoritmo durante a otimização, e o segundo com o valor da função custo a cada iteração.
plot_results(f_esfera, xi, yi, method = '1')
A Função Booth é uma função de $2$ dimensões definida como:
$$ $$\begin{equation} f(\mathbf{x}) = (x_{1} + 2x_{2} - 7)^{2} + (2x_{1} + x_{2} - 5)^{2} \end{equation}$$ $$Primeiro, temos seu gráfico:
f_booth = lambda xk: (xk[0] + 2*xk[1] - 7)**2 + (2*xk[0] + xk[1] - 5)**2
f_booth_grad = nd.Gradient(f_booth)
f_booth_hessian = nd.Hessian(f_booth)
plot_function(f_booth, title = 'Função Booth', ndim = 3)
chute = [5, 5]
xi, yi, grad_norm = GradientDescent(f_booth, f_booth_grad, chute = chute, tol=1e-10)
df2_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
Chute inicial: y = 164, x = [5, 5]
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 61 ***** Primeira iteração ***** f([5 5]) = 164.0 ||∇f([5 5])|| = 76.4198927 ***** Última iteração: 61° ***** f([1. 3.]) = 0.0 ||∇f([1. 3.])|| = 0.0
Que atingiu o critério de parada após $61$ iterações do Algoritmo. Abaixo temos os resultados em gráficos:
plot_results(f_booth, xi, yi, dim = [np.linspace(-10, 15, 500), np.linspace(-10, 15, 500)], method = '1')
Neste outro caso temos a Função Styblinski-Tang, definida como:
$$ $$\begin{equation} f(\mathbf{x}) = \dfrac{1}{2} \sum_{i = 1}^{n}(x_{i}^{4}-16x_{i}^{2}+5x_{i}) \end{equation}$$ $$Que foi considerada com $n = 2$, ou seja, duas dimensões e avaliada em $x_{0} = [-3, -4]$
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
chute = [-3, -4]
xi, yi, grad_norm = GradientDescent(f_tang, f_tang_grad, chute = chute, tol=1e-7)
df3_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
df3_1
Chute inicial: y = -49.0, x = [-3, -4]
| x1 | x2 | f(x1,x2) | ||∇f|| | |
|---|---|---|---|---|
| 0 | -3.00000000 | -4.00000000 | -49.00000000 | 61.59951299 |
| 1 | -2.89062500 | -2.07812500 | -69.58226913 | 17.80632896 |
| 2 | -2.91834593 | -3.19067526 | -76.76195048 | 11.42561239 |
| 3 | -2.90221876 | -2.83398927 | -78.25061301 | 2.32193091 |
| 4 | -2.90363926 | -2.90653571 | -78.33217526 | 0.10402795 |
| 5 | -2.90352553 | -2.90328682 | -78.33233035 | 0.00855307 |
| 6 | -2.90353471 | -2.90355395 | -78.33233140 | 0.00068935 |
| 7 | -2.90353397 | -2.90353242 | -78.33233141 | 0.00005565 |
| 8 | -2.90353403 | -2.90353416 | -78.33233141 | 0.00000449 |
| 9 | -2.90353403 | -2.90353402 | -78.33233141 | 0.00000036 |
| 10 | -2.90353403 | -2.90353404 | -78.33233141 | 0.00000042 |
| 11 | -2.90353403 | -2.90353399 | -78.33233141 | 0.00000140 |
| 12 | -2.90353403 | -2.90353403 | -78.33233141 | 0.00000011 |
| 13 | -2.90353403 | -2.90353402 | -78.33233141 | 0.00000013 |
| 14 | -2.90353403 | -2.90353403 | -78.33233141 | 0.00000015 |
| 15 | -2.90353403 | -2.90353402 | -78.33233141 | 0.00000018 |
| 16 | -2.90353403 | -2.90353403 | -78.33233141 | 0.00000008 |
Que atingiu o critério de tolerância $\leq 10^{-7}$ após $11$ iterações do Algoritmo.
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 17 ***** Primeira iteração ***** f([-3 -4]) = -49.0 ||∇f([-3 -4])|| = 61.59951299 ***** Última iteração: 17° ***** f([-2.90353403 -2.90353403]) = -78.33233141 ||∇f([-2.90353403 -2.90353403])|| = 8e-08
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '1')
Outro exemplo é a Função Rosenbrock, que é uma função $n-dimensional$ unimodal (com apenas 1 mínimo) sendo que esse mínimo se encontra em um vale parabólico estreito $[2]$.
$$ $$A função é definida como:
$$ $$\begin{equation} f(\mathbf{x}) = \sum_{i=1}^{n-1} \big(100 (x_{i+1} - x_ {i}^{2})^{2} + (x_{i}-1)^{2}\big) \end{equation}$$ $$Foi considerada a função com $n = 2$ e a constante (100) que multiplica o primeiro termo da Função Rosenbrock foi removido para facilitar a visualização com valores menores.
$$ $$Os resultados abaixo referem-se ao Método do Gradiente calculado com chute inicial $x_{0} = [-5, -25]$.
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
chute = [-5, -25]
xi, yi, grad_norm = GradientDescent(f_rosenbrock, f_rosenbrock_grad, chute = chute, tol=1e-10)
df4_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 2536, x = [-5, -25] *************** Resultados *************** Número de iterações: 335 ***** Primeira iteração ***** f([ -5 -25]) = 2536.0 ||∇f([ -5 -25])|| = 1016.9287094 ***** Última iteração: 335° ***** f([1. 1.]) = 0.0 ||∇f([1. 1.])|| = 0.0
Atingindo o critério desejado após 335 iterações.
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 4, 1000), np.linspace(-200, 250, 1000)], method = '1')
Podemos perceber que o método oscilou muito até que fosse atingido o critério de parada pré-estabelecido. Variando o critério de tolerância para uma precisão muito maior, pode-se perceber que o método não converge, ou seja, essas oscilações permanecem e não atingem o desejado. Esse resultado será importante para se comparar com os outros Algoritmos como o Método de Newton e o BFGS, que serão apresentados abaixo e que são teoricamente melhores que o Método de Descida do Gradiente.
A Função McCormick, por sua vez, é uma função de 2 variáveis reais e também unimodal, definida como:
$$ $$\begin{equation} f(\mathbf{x}) = \sin{(x_{1}+x_{2})} + (x_{1}-x_{2})^{2} - 1.5x_{1} + 2.5x_{2} + 1 \end{equation}$$ $$Essa função foi avaliada com $x_{0} = [-7, 6]$.
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
chute = [-7, 6]
xi, yi, grad_norm = GradientDescent(f_cormick, f_cormick_grad, chute = [-7, 6], tol=1e-10)
df5_1 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 194.6585290151921, x = [-7, 6] *************** Resultados *************** Número de iterações: 45 ***** Primeira iteração ***** f([-7 6]) = 194.65852902 ||∇f([-7 6])|| = 39.62530073 ***** Última iteração: 45° ***** f([-0.54719755 -1.54719755]) = -1.91322295 ||∇f([-0.54719755 -1.54719755])|| = 0.0
Que satisfez o critério de parada com $\lVert \nabla f(x_{k}) \rVert \leq 10^{-10}$ após $46$ iterações.
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '1')
O Método de Newton, diferentemente do método apresentado acima, consiste em se obter uma direção $d_k$ tal que:
\begin{equation} d_k = - (\nabla^{2}f(\vec{x}))^{-1} \cdot \nabla f(\vec{x}) \end{equation}$$ $$Em que $\nabla^{2}f(\vec{x}) = \mathbb{H}(x)_{(ij)} \equiv \frac{\partial^{2} f}{\partial x_{i} \partial x_{j} }$ é a matriz Hessiana.
Um fato a se destacar é que longe do ponto, $\nabla^{2}f(x_{k})$ pode ser não definida positiva. Caso isso ocorra, a direção do Método do Gradiente deve ser empregado como direção de descida.
A partir de então, calcular o novo ponto $x_{k+1}$ em que a norma do gradiente $\lVert \nabla f(x_{k+1}) \rVert$ seja menor e o custo da função $f(x_{k+1})$ diminua. O algoritmo esquematizado está abaixo:
Algoritmo 3 - Método de Newton com condição de Armijo
Dado um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$), para $\nabla f(x_{k}) \neq 0$ e $\nabla^{2}f(x_{k})$ matriz não singular, realizar os seguintes passos:
1. Calcular a direção de descida $d_{k}$ resolvendo o sistema $\nabla^{2}f(x_{k}) \, d_{k} = - \nabla f(x_{k})$;
2. Verificar se $d_{k}^{T}\nabla f(x_{k}) \leq -0.001 \lVert d_{k} \rVert \lVert \nabla f(x_{k}) \rVert$. Se não for satisfeita, usar $d_{k} = - \nabla f(x_{k})$ (Método do Gradiente);
3. Realizar Busca de Armijo e encontrar um $\lambda_k$ que satisfaça a Desigualdade 1;
4. Calcular o novo ponto $x_{k+1} = x_{k} + \lambda_{k} d_{k}$;
Em que tais etapas são realizadas com os mesmos critérios de parada dados no Método de Descida do Gradiente.
Ademais, a convergência do Método de Newton é quadrática, ou seja, a sequência {$x_k$} converge para $x^{*}$ quando:
\begin{equation} \lVert x_{k+1} - x^{*} \rVert \leq \epsilon \, \lVert x_{k} - x^{*} \rVert^{2} \end{equation}$\forall k \in \mathbb{N}$ e com $\epsilon > 0$.
Sendo assim teoricamente melhor que o método acima e o método apresentado a seguir.
A chamada para a Função que implementa o Método de Newton com Condição de Armijo é:
Newton(f, f_grad, f_hessian, chute, sigma, tol)
Onde:
$$ $$f é a função desejada (tipo: callable)
f_grad é o gradiente da função f (tipo: callable)
f_hessian é a hessiana da função f (tipo: callable)
chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)
sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)
tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = 1e-8 = $10^{-8}$)
E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.
$$ $$x são pontos $x_k$ calculados pelo Método
y são valores de $f(x_{k})$ calculados pelo Método
grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método
Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:
$$ $$x, y, _ = Newton(f, f_grad, f_hessian, chute, sigma, tol)
Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 3:
def Newton(f, f_grad, f_hessian, chute, sigma = 0.02, tol=1e-8):
"""Algoritmo de Descida do Gradiente com busca de Armijo. Esquema:
I - Direção de descida: d_k := −∇f(x).
II - Determinação do passo com busca de Armijo.
III - Obtém o próximo candidato.
Parâmetros
--------------------
f : callable
Função custo
f_grad : callable
Gradiente da função f
f_hessian : callable
Hessiana da função f
chute : array
Valor inicial de x ("chute")
sigma : float, opcional
Constante de Decréscimo de Armijo. Padrão = 0.02
tol : float, opcional
Tolerância desejada. Padrão = 1e-8
Saída
--------------------
xk : array
valores de xk
yk : array
valores da função custo e cada xk
grad_f : array
valores da norma do gradiente da f em xk"""
# Valores iniciais de xk, fk, grad_fk e max_iter
xk = chute
fk = f(xk)
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
max_iter = 500
# Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
num_iter = 0
x_pontos = [xk]
y_pontos = [fk]
grad_f = [grad_fk_norm]
f_hessian_xk = f_hessian(xk)
# Calcula nova iteração com busca de Armijo
while (grad_fk_norm > tol and num_iter < max_iter):
f_hessian_xk = f_hessian(xk)
if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16) or (np.linalg.norm(f_hessian_xk) > 1e16):
print("\nErro: Overflow\n")
break
# Não invertível
if abs(np.linalg.det(f_hessian_xk)) <= 1e-3:
break
# Determina a direção resolvendo o sistema: \nabla^{2} f(x_{k}) . p_{k} = - \nabla f(x_{k})
pk = -np.linalg.solve(f_hessian_xk, grad_fk)
lambda_ = 1
# Se estiver longe do ponto, a matriz Hessiana pode ser não definida positiva.
# Se isso acontecer, seguir pelo Método do Gradiente:
if not (np.dot(pk, grad_fk) < (-0.001 * np.linalg.norm(f_grad(xk)) * np.linalg.norm(pk))):
# Descida do gradiente
pk = -f_grad(xk)
# Faz a busca de Armijo, obtendo o passo lambda e a função custo naquele passo
lambda_, fk = Armijo_Search(f, xk, pk, sigma=sigma)
else:
lambda_, fk = Armijo_Search(f, xk, pk, sigma=sigma)
# calcula x_{k+1}
xk = xk + lambda_ * pk
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
# Itera mais uma vez
num_iter += 1
x_pontos.append(xk)
y_pontos.append(f(xk))
grad_f.append(grad_fk_norm)
# número máximo de iterações
if num_iter == max_iter:
print('\nNúmero máximo de iterações atingido.\n')
break
if abs(np.linalg.det(f_hessian_xk)) <= 1e-3:
print("\nMatriz Hessiana não invertível. Mudar o ponto inicial.\n")
return [], [], []
elif np.min(np.linalg.eigvals(f_hessian_xk)) <= 1e-3:
print("\nMatriz Hessiana não PD. Mudar o ponto inicial. Ponto de sela encontrado.\n")
return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)
Nesta parte algumas funções exemplo serão testadas utilizando o Método de Newton com busca linear de Armijo, como os exemplos 3, 4 e 5 do Método de Descida do Gradiente.
Testando o Método de Newton na Função Esfera definida anteriormente, temos:
f_esfera = lambda xk: xk[0]**2 + xk[1]**2
f_esfera_grad = nd.Gradient(f_esfera)
f_esfera_hessian = nd.Hessian(f_esfera)
plot_function(f_esfera, title = 'Função Esfera', ndim = 3)
chute = [4, -3]
xi, yi, grad_norm = Newton(f_esfera, f_esfera_grad, f_esfera_hessian, chute = chute, tol=1e-10)
df1_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
df1_2
| x1 | x2 | f(x1,x2) | ||∇f|| | |
|---|---|---|---|---|
| 0 | 4.00000000 | -3.00000000 | 25.00000000 | 10.00000000 |
| 1 | -0.00000000 | 0.00000000 | 0.00000000 | 0.00000000 |
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 2 ***** Primeira iteração ***** f([ 4 -3]) = 25.0 ||∇f([ 4 -3])|| = 10.0 ***** Última iteração: 2° ***** f([-0. 0.]) = 0.0 ||∇f([-0. 0.])|| = 0.0
O Método de Newton convergiu após 2 iterações. Abaixo temos os $\textit{plots}$ com mais alguns resultados:
plot_results(f_esfera, xi, yi, method = '2')
Para a Função Booth temos os seguintes resultados:
f_booth = lambda xk: (xk[0] + 2*xk[1] - 7)**2 + (2*xk[0] + xk[1] - 5)**2
f_booth_grad = nd.Gradient(f_booth)
f_booth_hessian = nd.Hessian(f_booth)
plot_function(f_booth, title = 'Função Booth', ndim = 3)
chute = [5, 5]
xi, yi, grad_norm = Newton(f_booth, f_booth_grad, f_booth_hessian, chute = chute, tol=1e-10)
df2_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
df2_2
| x1 | x2 | f(x1,x2) | ||∇f|| | |
|---|---|---|---|---|
| 0 | 5.00000000 | 5.00000000 | 164.00000000 | 76.41989270 |
| 1 | 1.00000000 | 3.00000000 | 0.00000000 | 0.00000000 |
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 2 ***** Primeira iteração ***** f([5 5]) = 164.0 ||∇f([5 5])|| = 76.4198927 ***** Última iteração: 2° ***** f([1. 3.]) = 0.0 ||∇f([1. 3.])|| = 0.0
Que atingiu o critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$ após 2 iterações. Aqui já vemos uma diferença considerável se comparada com o Método de Descida do Gradiente, uma vez que a convergência só foi atingida após $61$ iterações, enquanto que no Método de Newton foi praticamente imediata sob as mesmas condições de partida e parada.
plot_results(f_booth, xi, yi, dim = [np.linspace(-10, 15, 500), np.linspace(-10, 15, 500)], method = '2')
Neste caso temos novamente a Função Styblinski-Tang com $n = 2$ com um caso bastante interessante do Método de Newton. Primeiro, o algoritmo será executado com um passo inicial (chute inicial) de $x_{0} = [0, 2]$.
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
chute = [0, 2]
xi, yi, grad_norm = Newton(f_tang, f_tang_grad, f_tang_hessian, chute = chute, tol=1e-7)
df3_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Matriz Hessiana não PD. Mudar o ponto inicial. Ponto de sela encontrado. *************** Resultados *************** Número de iterações: 5 ***** Primeira iteração ***** f([0 2]) = -19.0 ||∇f([0 2])|| = 13.72953022 ***** Última iteração: 5° ***** f([0.15673126 2.74680277]) = -24.8338343 ||∇f([0.15673126 2.74680277])|| = 0.0
O Algoritmo 3 identificou que a condição sobre a Matriz Hessiana ser Positiva Definida foi violada. Trata-se de um ponto de sela.
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '2')
Esse é um clássico problema do Método de Newton, tanto com $1$ variável real como multivariado, a escolha incial é importante para a convergência do método.
$$ $$Mudando o ponto de começo $x_{0}$ para $x_{0} = [-3, -4]$, temos os seguintes resultados:
chute = [-3, -4]
xi, yi, grad_norm = Newton(f_tang, f_tang_grad, f_tang_hessian, chute = [-3, -4], tol=1e-7)
df3_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 6 ***** Primeira iteração ***** f([-3 -4]) = -49.0 ||∇f([-3 -4])|| = 61.59951299 ***** Última iteração: 6° ***** f([-2.90353403 -2.90353403]) = -78.33233141 ||∇f([-2.90353403 -2.90353403])|| = 0.0
Chegando no mesmo valor que o Método do Gradiente com apenas 6 passos.
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '2')
Que convergiu para o mínimo da função.
A função testada é novamente a Função Rosenbrock com n = 2 utilizada no Exemplo do Método anterior.
Os resultados abaixo referem-se ao Método de Newton calculado com chute inicial $x_{0} = [-5, -25]$.
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
f_rosenbrock_hessian = nd.Hessian(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
chute = [-5, -25]
xi, yi, grad_norm = Newton(f_rosenbrock, f_rosenbrock_grad, f_rosenbrock_hessian, chute = chute, tol=1e-10)
df4_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 13 ***** Primeira iteração ***** f([ -5 -25]) = 2536.0 ||∇f([ -5 -25])|| = 1016.9287094 ***** Última iteração: 13° ***** f([1. 1.]) = 0.0 ||∇f([1. 1.])|| = 0.0
No Método do Gradiente, para esse problema, o critério de parada $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$ foi satisfeito somente após $335$ iterações. O Método de Newton, por sua vez, sob as mesmas condições, atingiu o critério após apenas $13$ iterações, sendo outra diferença notável da convergência dos dois métodos.
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 4, 1000), np.linspace(-200, 250, 1000)], method = '2')
Além disso, o Método de Newton não oscilou como o Método do Gradiente para encontrar a solução em que o critério fosse satisfeito, que acaba sendo a solução da Função Rosenbrock bi-dimensional.
A Função McCormick foi avaliada com $x_{0} = [-7, 6]$ e com critério $\lVert\nabla f(\mathbf{x}) \rVert \leq 10^{-10}$, como anteriormente.
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
f_cormick_hessian = nd.Hessian(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
chute = [-7, 6]
xi, yi, grad_norm = Newton(f_cormick, f_cormick_grad, f_cormick_hessian, chute = chute, tol=1e-10)
df5_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 5 ***** Primeira iteração ***** f([-7 6]) = 194.65852902 ||∇f([-7 6])|| = 39.62530073 ***** Última iteração: 5° ***** f([-0.54719755 -1.54719755]) = -1.91322295 ||∇f([-0.54719755 -1.54719755])|| = 0.0
Novamente encontrando os mesmos valores do Método de Newton em muito menos iterações do Algoritmo.
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '2')
A Função Colville é uma função de $4$ dimensões, definida como:
$$ $$\begin{equation} f(\mathbf x) = 100(x_{1}^2-x_{2})^2 + (x_{1}-1)^2 + (x_{3} -1)^2 + 90(x_{3}^2-x_{4})^2\\+10.1((x_{2} - 1)^2+ (x_{4}-1)^2) +19.8(x_{2} -1)(x_{4} -1) \end{equation}$$ $$Que será avaliada no ponto $x_{0} = [3, 1, 5, 4]$:
def f_colville(xk):
"""
Função Colville - 4 dimensões: https://www.sfu.ca/~ssurjano/colville.html
Parâmetros
--------------------
xk : np.array or list
array de pontos xk
Saída
--------------------
y : array
valores da função custo
"""
x1 = xk[0]
x2 = xk[1]
x3 = xk[2]
x4 = xk[3]
term1 = 100 * (x1**2-x2)**2;
term2 = (x1-1)**2;
term3 = (x3-1)**2;
term4 = 90 * (x3**2-x4)**2;
term5 = 10.1 * ((x2-1)**2 + (x4-1)**2);
term6 = 19.8*(x2-1)*(x4-1);
y = term1 + term2 + term3 + term4 + term5 + term6
return y
f_colville_grad = nd.Gradient(f_colville)
f_colville_hessian = nd.Hessian(f_colville)
chute = [3, 1, 5, 4]
xi, yi, grad_norm = Newton(f_colville, f_colville_grad, f_colville_hessian, chute = chute, tol=1e-8)
df6_2 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'x3':xi[:,2], 'x4':xi[:,3], 'f(x1,x2,x3,x4)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'x3', 'x4', 'f(x1,x2,x3,x4)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
*************** Resultados *************** Número de iterações: 15 ***** Primeira iteração ***** f([3 1 5 4]) = 46200.9 ||∇f([3 1 5 4])|| = 39215.92871169 ***** Última iteração: 15° ***** f([1. 1. 1. 1.]) = 0.0 ||∇f([1. 1. 1. 1.])|| = 0.0
df6_2
| x1 | x2 | x3 | x4 | f(x1,x2,x3,x4) | ||∇f|| | |
|---|---|---|---|---|---|---|
| 0 | 3.00000000 | 1.00000000 | 5.00000000 | 4.00000000 | 46200.90000000 | 39215.92871169 |
| 1 | 2.29966013 | 2.93255387 | 4.49767692 | 17.87089851 | 4627.45416071 | 4399.11465877 |
| 2 | 1.66303161 | 1.70962687 | 3.09191394 | 6.84928366 | 1210.45916989 | 3123.64881103 |
| 3 | 1.02829368 | 0.25151366 | 2.56382101 | 5.83413244 | 286.58824030 | 767.13016041 |
| 4 | 0.82740311 | 0.48596161 | 1.92643823 | 3.12318333 | 62.53754579 | 421.28883716 |
| 5 | 0.54946097 | 0.15520747 | 1.64994513 | 2.56335455 | 10.79500134 | 102.82494463 |
| 6 | 0.50791572 | 0.24068237 | 1.41680691 | 1.93193006 | 1.54159644 | 40.57739563 |
| 7 | 0.58543341 | 0.33528193 | 1.30739964 | 1.69270586 | 0.48876832 | 8.82368418 |
| 8 | 0.79961196 | 0.59452106 | 1.19068036 | 1.40342568 | 0.36160394 | 18.08480041 |
| 9 | 0.83862492 | 0.70294929 | 1.13752332 | 1.29113716 | 0.08062932 | 1.56255741 |
| 10 | 0.99250921 | 0.96141233 | 1.02560337 | 1.03918358 | 0.07177742 | 11.74202871 |
| 11 | 0.99147987 | 0.98296381 | 1.00919960 | 1.01806179 | 0.00030743 | 0.17871781 |
| 12 | 0.99982048 | 0.99957110 | 1.00026150 | 1.00044089 | 0.00000127 | 0.04548513 |
| 13 | 0.99999827 | 0.99999652 | 1.00000180 | 1.00000352 | 0.00000000 | 0.00003561 |
| 14 | 1.00000000 | 1.00000000 | 1.00000000 | 1.00000000 | 0.00000000 | 0.00000000 |
Convergindo para o mínimo global $f(1,1,1,1) = 0$ após $15$ passos.
O método $\textbf{B}$royden–$\textbf{F}$letcher–$\textbf{G}$oldfarb–$\textbf{S}$hanno ($\textbf{BFGS}$) é um método Quasi-Newton, sendo uma alternativa apresentada para o problema encontrado no Método de Newton apresentado acima, em que é necessário calcular a matriz Hessiana $\nabla^{2}{f(\vec{x})}$, ($x \in \Re^{n}$).
Tal problema é contornado pelo cálculo de uma matriz que serve como uma aproximação para a matriz Hessiana, não sendo necessários, portanto, o seu cálculo diretamente e o cálculo de sua inversa, como foi realizado no método acima.
O algoritmo e sua esquematização foram inspiradas em $[4]$.
Algoritmo 4 - Método BFGS com condição de Armijo
A partir de um ponto $x_{0} \in \Re^{n}$ ($\textit{initial guess}$) e uma matriz inicial $B_{0} \equiv \mathbb{I}_{n}$, com $n \in \mathbb{N}$, realizar os seguintes passos:
$$ $$1. Obter uma direção de descida $d_{k} = - B_{k} \cdot \nabla f(x_{k})$;
$$ $$2. Realizar busca linear e determinar $\lambda_{k}$ que satisfaz a condição de Armijo (Desigualdade 1);
$$ $$3. Fazer $s_{k} = \lambda_{k} \, d_{k}$ e atualizar o ponto $x_{k+1} = x_{k} + \lambda_{k} \, d_{k} = x_{k} + s_{k}$;
$$ $$4. Calcular a diferença de gradientes $y_{k} = \nabla f(x_{k+1}) - \nabla f(x_{k})$;
$$ $$5. Atualizar a matriz $B_{k+1}$ que aproxima a Hessiana:
\begin{equation}\tag{2} B_{k+1} = B_{k} + \dfrac{(s^{T}_{k} y_{k}+ y^{T}_{k}B_{k}y_{k})(s_{k}s^{T}_{k})}{(s^{T}_{k}y_{k})^{2}} - \dfrac{B_{k}y_{k}s^{T}_{k} + s_{k}y^{T}_{k}B_{k}}{s^{T}_{k}y_{k}} \end{equation}$$ $$Em que tais etapas são realizadas até os critérios apresentados no começo deste trabalho forem satisfeitos.
$$ $$Vale ressaltar ainda que a convergência do Método BFGS é superlinear, ou seja, a sequência {$x_n$} converge para $x^{*}$ quando:
\begin{equation} \lim_{n→\infty} \dfrac{\lVert x_{n+1} - x^{*} \rVert}{\lVert x_{n} - x^{*} \rVert} = 0 \end{equation}$$ $$Sendo assim, teoricamente melhor que o Método de Descida do Gradiente, porém pior que o Método de Newton.
A chamada para a Função que implementa o Método BFGS com Condição de Armijo é:
BFGS(f, f_grad, chute, sigma, tol)
Onde:
$$ $$f é a função desejada (tipo: callable)
f_grad é o gradiente da função f (tipo: callable)
chute é o $x_{0}$, ponto inicial (tipo: np.array ou list)
sigma é a constante de decréscimo de Armijo (tipo: float $\in (0, 1)$. Padrão $= 0.02$)
tol é a tolerância desejada ($\lVert \nabla f(x_{k+1}) \rVert \leq 10^{-8}$) por exemplo (tipo: float $> 0$. Padrão = $10^{-8}$)
E pode retornar até 3 argumentos que são tipo np.array utilizados para fazer gráficos e outras análises visuais.
$$ $$x são pontos $x_k$ calculados pelo Método
y são valores de $f(x_{k})$ calculados pelo Método
grad_f são valores de $\lVert \nabla f(x_k) \rVert$ calculados pelo Método
Para suprimir algum retorno, por exemplo, o grad_f, basta fazer:
$$ $$x, y, _ = BFGS(f, f_grad, chute, sigma, tol)
Exemplos de como utilizar tais funções serão fornecidos mais adiante com funções pré-estabelecidas. Abaixo temos a implementação do Algoritmo 4:
def BFGS(f, f_grad, chute, sigma = 0.02, tol=1e-8):
""" Algoritmo Broyden–Fletcher–Goldfarb–Shanno (BFGS) com busca de Armijo.
Parâmetros
--------------------
f : callable
Função custo
f_grad : callable
Gradiente da função f
chute : array
Valor inicial de x ("chute")
sigma : float, opcional
Constante de decréscimo de Armijo. Padrão = 0.02
tol : float, opcional
Tolerância desejada. Padrão = 1e-8
Saída
--------------------
xk : array
valores de xk
yk : array
valores da função custo em cada xk
grad_f : array
valores da norma do gradiente da f em xk"""
# dimensão do problema e matriz B_{k = 0} = Id(n) (identidade de ordem n)
n = len(chute)
Bk = np.eye(len(chute))
max_iter = 500
# Valores iniciais de x_{k+1}, x_{k}, f_{x_{k}}, grad_{f(x_{k})} e norma do gradiente
xk = chute
xprev = chute
fk = f(xk)
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
# Inicializa o número de iterações e a lista para fazer os plots dos valores de x e y
num_iter = 0
x_pontos = [xk]
y_pontos = [fk]
grad_f = [grad_fk_norm]
print(f'Chute inicial: y = {fk}, x = {xk} \n')
while (grad_fk_norm > tol and num_iter < max_iter):
if (np.linalg.norm(f(xk)) > 1e16) or (grad_fk_norm > 1e16):
print("\nErro: Overflow\n")
break
# direção (método de Newton)
dk = - Bk @ f_grad(xprev)
# performa busca de Armijo com x_{k} (x prévio)
lambda_, _ = Armijo_Search(f, xprev, dk, sigma)
# direção e valor de x_{k+1} = x_{k} + \lambda * d_{k} = x_{k} + s_{k}
sk = lambda_ * dk
xk = xk + lambda_ * dk
# y_{k} é a diferença do gradiente da f em x_{k+1} com o gradiente da f em x_{k}
yk = np.subtract(f_grad(xk), f_grad(xprev))
# arruma os tamanhos
yk = np.array([yk]).reshape(-1, 1)
sk = np.array([sk]).reshape(-1, 1)
# denominador do BFGS update
denominador = np.dot(sk.T, yk)[0] # escalar
# numerador da segunda parcela do BFGS update
num_1 = ((denominador + (yk.T @ (Bk @ yk))) * (sk @ sk.T))/(denominador ** 2)
# numerador da terceira parcela do BFGS update
num_2 = ((Bk @ (yk @ sk.T)) + (sk @ (yk.T @ Bk)))/denominador
# BFGS update
Bk = Bk + np.subtract(num_1, num_2)
grad_fk = f_grad(xk)
grad_fk_norm = np.linalg.norm(grad_fk)
x_pontos.append(xk)
y_pontos.append(np.round(f(xk), 10))
grad_f.append(grad_fk_norm)
xprev = xk
num_iter += 1
# número máximo de iterações
if num_iter == max_iter:
print('\nNúmero máximo de iterações atingido.\n')
return np.array(x_pontos), np.array(y_pontos), np.array(grad_f)
Nesta parte algumas funções exemplo serão testadas utilizando o Método BFGS com busca linear de Armijo.
f_tang = lambda xk: (xk[0]**4 - 16*xk[0]**2 + 5*xk[0] + xk[1]**4 - 16*xk[1]**2 + 5*xk[1])/2
f_tang_grad = nd.Gradient(f_tang)
f_tang_hessian = nd.Hessian(f_tang)
plot_function(f_tang, title = 'Função Styblinski-Tang com n = 2', ndim = 3, angle = (40, 50))
chute = [-3, -2]
xi, yi, grad_norm = BFGS(f_tang, f_tang_grad, chute = chute, tol=1e-7)
df1_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = -68.0, x = [-3, -2] *************** Resultados *************** Número de iterações: 11 ***** Primeira iteração ***** f([-3 -2]) = -68.0 ||∇f([-3 -2])|| = 18.82817038 ***** Última iteração: 11° ***** f([-2.90353403 -2.90353403]) = -78.33233141 ||∇f([-2.90353403 -2.90353403])|| = 6e-08
plot_results(f_tang, xi, yi, dim = [np.linspace(-5, 5, 500), np.linspace(-5, 5, 500)], method = '3')
Nessa primeira função é interessante perceber que o número de iterações para o critério de parada foi maior que no Método de Newton.
Testando a função Rosenbrock $2$-$d$ no método BFGS, temos:
f_rosenbrock = lambda xk: (xk[1] - xk[0]**2)**2 + (1 - xk[0])**2
f_rosenbrock_grad = nd.Gradient(f_rosenbrock)
f_rosenbrock_hessian = nd.Hessian(f_rosenbrock)
plot_function(f_rosenbrock, title = 'Função Rosenbrock $n=2$', ndim = 3, dom = np.linspace(-2, 2, 500), angle = (30, 40))
chute = [-5, -25]
xi, yi, grad_norm = BFGS(f_rosenbrock, f_rosenbrock_grad, chute = [-5, -25], tol=1e-10)
df2_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 2536, x = [-5, -25] *************** Resultados *************** Número de iterações: 17 ***** Primeira iteração ***** f([ -5 -25]) = 2536.0 ||∇f([ -5 -25])|| = 1016.9287094 ***** Última iteração: 17° ***** f([1. 1.]) = 0.0 ||∇f([1. 1.])|| = 0.0
Percebe-se que o Método BFGS teve mais iterações que o Método de Newton e menos iterações que o Método do Gradiente, até que o critério de parada fosse atingido.
plot_results(f_rosenbrock, xi, yi, dim = [np.linspace(-13, 15, 1000), np.linspace(-200, 250, 1000)], method = '3')
f_cormick = lambda xk: np.sin(xk[0] + xk[1]) + (xk[0] - xk[1])**2 - 1.5*xk[0] + 2.5*xk[1] + 1
f_cormick_grad = nd.Gradient(f_cormick)
f_cormick_hessian = nd.Hessian(f_cormick)
plot_function(f_cormick, title = 'Função McCormick', ndim = 3, dom = np.linspace(-4, 4, 500), angle = (30, 50))
chute = [-7, 6]
xi, yi, grad_norm = BFGS(f_cormick, f_cormick_grad, chute = chute, tol=1e-10)
df3_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'f(x1,x2)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'f(x1,x2)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 194.6585290151921, x = [-7, 6] *************** Resultados *************** Número de iterações: 9 ***** Primeira iteração ***** f([-7 6]) = 194.65852902 ||∇f([-7 6])|| = 39.62530073 ***** Última iteração: 9° ***** f([-0.54719755 -1.54719755]) = -1.91322296 ||∇f([-0.54719755 -1.54719755])|| = 0.0
plot_results(f_cormick, xi, yi, dim = [np.linspace(-45, 8, 1000), np.linspace(-45, 8, 1000)], method = '3')
Pelo Método BFGS, o mínimo da função foi encontrado com apenas 8 iterações, o que não ocorreu no Método de Newton, que resultou em um valor diferente do esperado. Se comparado com o Método do Gradiente, o número de iterações foi muito menor (46 iterações foram necessárias no primeiro Método), o que teoricamente faz sentido uma vez que o tipo de convergência é melhor se comparado com o Método do Gradiente.
Testando agora na Função Colville:
f_colville_grad = nd.Gradient(f_colville)
chute = [3, 1, 5, 4]
xi, yi, grad_norm = BFGS(f_colville, f_colville_grad, chute = chute, tol=1e-8)
df4_3 = pd.DataFrame(data = {'x1':xi[:,0], 'x2':xi[:,1], 'x3':xi[:,2], 'x4':xi[:,3], 'f(x1,x2,x3,x4)':yi, '||∇f||':grad_norm}, index = pd.Index([i for i in range(len(xi))]), columns = ['x1', 'x2', 'x3', 'x4', 'f(x1,x2,x3,x4)', '||∇f||'])
printa_resultados(chute, xi, yi, grad_norm)
Chute inicial: y = 46200.9, x = [3, 1, 5, 4] *************** Resultados *************** Número de iterações: 77 ***** Primeira iteração ***** f([3 1 5 4]) = 46200.9 ||∇f([3 1 5 4])|| = 39215.92871169 ***** Última iteração: 77° ***** f([1. 1. 1. 1.]) = 0.0 ||∇f([1. 1. 1. 1.])|| = 0.0
df4_3
| x1 | x2 | x3 | x4 | f(x1,x2,x3,x4) | ||∇f|| | |
|---|---|---|---|---|---|---|
| 0 | 3.00000000 | 1.00000000 | 5.00000000 | 4.00000000 | 46200.90000000 | 39215.92871169 |
| 1 | 0.65527344 | 1.37612305 | -4.23046875 | 4.90805664 | 15485.72247885 | 19923.11661312 |
| 2 | -2.34055222 | 1.70232693 | -3.63569501 | 7.61445932 | 4823.45883370 | 8221.27443435 |
| 3 | -1.23940941 | 3.33747143 | -3.66775263 | 11.54402139 | 2345.12224793 | 2751.91227485 |
| 4 | -0.40102996 | 0.00289211 | -3.86812964 | 13.23333383 | 1577.26670286 | 2426.65070686 |
| ... | ... | ... | ... | ... | ... | ... |
| 72 | 0.99997852 | 0.99995828 | 1.00001805 | 1.00003609 | 0.00000000 | 0.00056475 |
| 73 | 0.99999924 | 0.99999858 | 1.00000052 | 1.00000093 | 0.00000000 | 0.00006695 |
| 74 | 1.00000000 | 1.00000001 | 0.99999999 | 0.99999998 | 0.00000000 | 0.00000135 |
| 75 | 1.00000000 | 1.00000000 | 1.00000000 | 1.00000000 | 0.00000000 | 0.00000001 |
| 76 | 1.00000000 | 1.00000000 | 1.00000000 | 1.00000000 | 0.00000000 | 0.00000000 |
77 rows × 6 columns
O mínimo da função foi encontrado com $78$ passos, mais que no Método de Newton.
Nessa parte de testes, o primeiro passo é escrever a função que deseja minimizar. Para isso, siga:
$$ $$f = lambda x: x[0]**2 + x[1]**2 + 2*x[0]*x[1] + x[2]
Ou
$$ $$def f(x):
return x[0]**2 + x[1]**2 + 2*x[0]*x[1] + x[2]
Para minimizar a função $f(x, y, z) = x^{2} + y^{2} + 2xy + z$, por exemplo.
$$ $$Depois escolha o método a ser utilizado com um número inteiro (1, 2 ou 3):
Método do Gradiente = 1
Método de Newton = 2
Método BFGS = 3
Finalmente, escolha o chute inicial (uma lista com a mesma dimensão da função), a constante de decréscimo de Armijo $\sigma$ (float $> 0$) e a tolerância (1e-8 para $10^{-8}$), por exemplo.
Função Objetivo:
$$ $$\begin{equation} f(\mathbf{x}) = (x_{1} + 2x_{2} - 7)^{2} + (2x_{1} + x_{2} - 5)^{2} + e^{-2y} \end{equation}$$ $$Entrada da forma, para o Método do Gradiente:
$$ $$f = lambda x: (x[0] + 2*x[1] - 7)**2 + (2*x[0] + x[1] - 5)**2 + np.exp(-2*x[1])
minimiza_f(f, metodo = 1, chute = [8, 15], sigma = 0.02, tol = 1e-10)
Cujo resultado está abaixo:
f = lambda x: (x[0] + 2*x[1] - 7)**2 + (2*x[0] + x[1] - 5)**2 + np.exp(-2*x[1])
metodo = 2
chute = [8, 15]
sigma = 0.5
tol = 1e-10
minimiza_f(f = f, metodo = metodo, chute = chute, sigma = sigma, tol = tol)
Método de Newton com condição de Armijo *************** Resultados *************** Número de iterações: 24 ***** Primeira iteração ***** f([ 8 15]) = 1637.0 ||∇f([ 8 15])|| = 241.93387526 ***** Última iteração: 24° ***** f([0.99890135 3.00137331]) = 0.00247535 ||∇f([0.99890135 3.00137331])|| = 0.0
[1] SURJANOVIC, S., BINGHAM, D., Simon Fraser University. Virtual Library of Simulation Experiments: Test Functions and Datasets. [https://www.sfu.ca/~ssurjano/optimization.html]
$$ $$[2] Test functions for optimization, Wikipedia. [https://en.wikipedia.org/wiki/Test_functions_for_optimization]
$$ $$[3] Real Valued Test Functions, HeuristicLab. [https://dev.heuristiclab.com/trac.fcgi/wiki/Documentation/Reference/Test%20Functions] $$ $$
[4] Broyden–Fletcher–Goldfarb–Shanno algorithm, Wikipedia.[https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm]
$$ $$[5] FRIEDLANDER, A. Elementos de Programação Não-Linear. Universidade Estadual de Campinas.[https://www.ime.unicamp.br/~friedlan/livro.pdf]